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Abstract — Ellipse and ellipsoid fitting has been extensively re- 
searched and widely applied. Although traditional fitting methods 
provide accurate estimation of ellipse parameters in the low- 
noise case, their performance is compromised when the noise 
level or the ellipse eccentricity are high. A series of robust fitting 
algorithms are proposed that perform well in high-noise, high- 
eccentricity ellipse/spheroid (a special class of ellipsoid) cases. 
The new algorithms are based on the geometric definition of an 
ellipse/spheroid, and improved using global statistical properties 
of the data. The efficacy of the new algorithms is demonstrated 
through simulations. 

I. Introduction 

Ellipse and ellipsoid fitting has been extensively studied and 
has broad applications. For example, an ellipse can serve as 
a geometric primitive model in computer vision and pattern 
recognition, which allows reduction and simplification of data 
for higher level processing. Ellipse fitting is also essential in 
applied sciences from observational astronomy (orbit estima- 
tion) to structural geology (strain in natural rocks). Moreover, 
ellipsoid fitting (e.g., minimum volume ellipsoid estimator 
[10]) is a useful tool in outlier detection. 

In this paper, we consider the following issues; In Section [ll] 
the ellipse fitting problem is formulated and several important 



algorithms are reviewed. In Section III a new objective func- 
tion based on the geometric definition of ellipse is introduced, 
and it is further extended to three ellipse fitting algorithms in 
Section IV] Then a spheroid fitting algorithm is proposed in 
Section W which is followed by the experimental results in 
Section El 

II. Formulation and Background 

As a classical problem, ellipse fitting has a rich literature. 
Various algorithms have been proposed from very different 
perspectives. We start our discussion by reviewing several 
classes of the most important ellipse fitting algorithms. 

The problem of ellipse fitting can be modeled as follows. 
We have data points {zi = (xi, y^)}"^]^, which are points on 
an ellipse corrupted by noise. The objective is to fit an ellipse 
with unknown parameters f3 to the data points, so that the total 
error is minimized. The measure of error differs for different 
classes of algorithms. 
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The most intuitive class of ellipse fitting algorithms is 
algebraic fitting, which uses algebraic error as a measure 
of error An ellipse can be described by an implicit function 
P(/9) = Ax'^ + Bxy + Ci/ + Dx + Ey + F = 0, where 
/3 — {A, B, C, D, E, F) denotes the ellipse parameters. The 
algebraic error from a data point Zi = {xi,yi) to the ellipse 
is thus Axf + Bxiyi + Cyf + Dx^ + Ey^ + F. 

The most efficient and widely used algorithm in this cate- 
gory was proposed by Fitzgibbon et al. [7]. It is a least squares 
optimization problem with respect to the algebraic criterion; 

n 

min ^(Ax,^ + Bxiyi + Cyf + Dx^ + Ey, + Ff{\) 

^ 1=1 



s.t. 



B^ -AAC = 1. 



The constraint ensures that the optimization result is an ellipse 
instead of a general conic section and also prevents the 
problem of parameter free scaling. This optimization problem 
is further reduced to a generalized eigenvalue problem, which 
can be efficiently solved. In addition, Halir et al. improved 
the algorithm into a numerically stable version [8]. Algebraic 
fitting has the advantage of low computational complexity. 
However, the algebraic criterion lacks geometric interpretation, 
and the algorithm is difficult to generalize to three dimensions 
due to its non-linear constraint. 

To overcome the shortcomings of algebraic fitting, Ahn et al. 
proposed orthogonal least squares fitting (OLSF) [2]. OLSF 
minimizes the sum of squared Euclidian distance, which is 
defined as the orthogonal distance from the data point to the 
ellipse; 



min > 
p ^ 

i=i 



(2) 



where z^{j3) is the orthogonal contacting point, which is 
the point on the ellipse that has the shortest distance to 
the corresponding data point. OLSF has a clear geometric 
interpretation and features high accuracy. Moreover, it can be 
generalized to the three-dimensional case [1]. Unfortunately, 
OLSF is computationally intensive. It employs the iterative 
Gauss-Newton algorithm, and in each iteration, the orthogonal 
contacting points have to be found for each data point, being 
iterative itself. 

Various extensions to OLSF algorithms have also been 
proposed. Angular information is incorporated into the OLSF 



algorithm in Watson's 2002 paper [12], in which the orthog- 
onal geometric distance is replaced by the distance along the 
known measurement angle. Moreover, instead of the I2 norm 
in OLSF, li, loo and Ip norms have been considered as well 
[11] [3] [4]. 

The third class of algorithms consists of Maximum Likeli- 
hood (ML) algorithms, which were proposed in [5] and [9]. 
The key steps of the two ML algorithms, the fundamental nu- 
merical scheme [5] and the heteroscedastic errors-in-variables 
scheme [9], have been proven to be equivalent in [6]. 

In [5], Chojnacki et al. assume that the data points are 
independent and have a multivariate normal distribution: 
Zi ^ A/^(z^,AzJ. The ML solution is then reduced to an 
optimization problem based on Mahalanobis distance: 



(z,;-zK/3))^A,7(z,-zK/3)). 



(3) 



The fundamental numerical scheme is implemented to solve a 
variational equation iteratively by solving an eigenvalue prob- 
lem at each iteration until it converges. The ML algorithms are 
accurate with moderate computational cost. However, when 
the noise is large or the eccentricity of the ellipse is large, the 
algorithm breaks down, because when any data point is close 
to the center of the estimated ellipse, one of the matrices in 
the algorithm has elements that tend to infinity. 

All three classes of algorithms described above have their 
advantages and perform well when the noise level is low. 
However, they share a common disadvantage of not being 
robust enough for highly noisy data. In this paper, we propose 
a series of algorithms that are resistant to large noise, and 
can be generalized to three-dimensions easily, with competitive 
accuracy and moderate computational cost. 

III. A New Geometric Objective Function 

Note that both OLSF and ML algorithms estimate the 
nuisance parameters z^ (point on the ellipse that generates the 
data point) in addition to the parameters of the ellipse. It is 
desirable to bypass these nuisance parameters and have a more 
intrinsic fitting method with a clear geometric interpretation. 

Recall the geometric definition of an ellipse. An ellipse is 
the locus of points such that the sum of the distances from that 
point to two other fixed points (called the foci of the ellipse) 
is constant. I.e., z is a point on the ellipse if and only if 
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(4) 



where j] • j] denotes the I2 norm, Ci and C2 are the two foci, 
and a is the length of the semi-major axis. 

Based on the geometric definition, the ellipse fitting problem 
can be naturally formulated as an optimization problem with 
a new geometric objective function: 
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(5) 



where n denotes the number of data points. 

The objective function in (|5]l has several advantages. First, it 
has a clear geometric interpretation: it is the expected squared 




Fig. 1. Derivation of the expected contribution of a noisy data point to our 
objective function. 



difference between the sum of the distances from the data 
points to the foci and the length of the major axis. Second, the 
parameters in the objective function, Ci, C2 and a, are intrinsic 
parameters of ellipses, which are translation and rotation 
invariant. The third advantage is that the objective function is 
ellipse specific, and thus no extra constraints are needed. As 
a result, the minimization problem can be readily solved by 
gradient descent algorithms. Lastly, and most importantly, the 
objective function possess one more property that contributes 
to the robustness of our algorithm, which is demonstrated by 
comparing with the objective function of OLSF. 

Assume that the data points are independent and have 
a multivariate normal distribution: Zi ^ J\f{z^, A^-)- The 
expected contribution of a noisy data point to the objective 
function of OLSF is approximately 



Ef 



(6) 



Note that this is homogeneous for all the data points. However, 
for our objective function, the expected contribution can be 
approximated as (see Fig. [T]| 

E[/zJ = E (Aj/(sin6'ji +sin6lj2) + Ax(cose,i -cose,J)^ 
= 2(7^ + 2cr^(sin Oil sin 6*^2 — cos On cos 6*^2) 
= 2a\l-cos{9a + 0,2)) 

= 2a2(l + cosC), (7) 

where Q is the angle ZciZiC2. This quantity is heterogeneous 
around the periphery of the ellipse: the objective function puts 
a large weight on the data points located at two ends of the 
major axis. As a result, our algorithm provides a highly robust 
and accurate major axis orientation estimate. 

However, the objective function in (|5]l has a global minimum 
at infinity. When the foci move away from each other toward 
infinity and the semi-major axis length tends to infinity, the 
value of the objective function approaches zero. So when the 
noise level is high, the local minimum (that we desire in this 
case) is smeared out by the noise, resulting in the estimated 
foci slipping away along the major axis. 

In order to take advantage of the robustness of ellipse axis 
orientation estimation while overcoming the problem of having 
a global minimum at infinity, we propose three modified 
algorithms in the next section, so as to increase the robustness 
and accuracy of the algorithm directly using our new objective 
function. 



IV. Three Modified Algorithms for 
Ellipse Fitting 



its major axis coinciding with the x-axis). The resulting data 
points are 



In this section, we introduce three modified algorithms 



based on the objective function proposed in Section III Each 
of the three algorithms overcomes the problem of a global 
minimum at infinity and achieves a robust and accurate result. 

A. Penalized Objective Function 

The most intuitive way to eliminate the global minimum at 
infinity is penalization. The penalty term should favor smaller 
semi-major lengths, so that the foci will not diverge. As a 
result, we propose the following penalized objective function: 



1 



|zi-Ci|| + ||zi-C2||-2a)^+AainaxO-exp 



i=l 




where the second term is the penalty term, A is a tuning 
parameter, Omax is an estimated upper bound for the semi- 
major length a, and a is an estimate of the noise standard 
deviation, a and Omax are estimated during the initialization 
procedure. 

The term exp((^^-S— )^) tends to infinity rapidly when a 
exceeds Omax, thus eliminating the global minimum at infinity. 
The penalty term is also proportional to the estimated noise 
level. When the noise level is high, the penalty term is large, 
so as to increase the robustness of the algorithm. On the other 
hand, when the noise level is low, the penalty term is small 
to ensure a small bias added to the objective function. The 
coefficient flmax is a scaling factor to accommodate the size 
of the ellipse we are fitting. 

In the initialization, a,„ax is estimated as the largest distance 
from a data point to the mean of all the data points Zmean- a is 
initialized as the mean of the distance from the data points to 
Zmean- To estimate (T, we run the gradient descent algorithm for 
Q and terminate when a exceeds Omax- Then a is estimated as 
the square root of the resulting objective function value. This 
simple noise level estimation method suffices for our purposes 
here. 

The penalized optimization problem (|8]l can be readily 
solved by gradient descent algorithms once initialized. The 
global minimum at infinity is eliminated and the resulting 
algorithm has good accuracy and high robustness, which we 
will show in Section [Vl] via simulation results. 

B. Axial Guided Ellipse Fitting 

The second algorithm we propose to overcome the problem 



described in Section III is axial guided fitting. Recall that the 
major axis orientation estimate is accurate and robust. So the 
only thing left to be determined is the size of the ellipse, i.e., 
the semi-major length a and semi-minor length b. 

To estimate a and b, we first solve ^ to estimate the major- 
axis orientation (j> and ellipse center {xc, j/c)- Then we translate 
and rotate the data points {{xi,yi)}"^^ so that the estimated 



and 



x[ = cos (t>{xi - Xc) + sin ^(yi - yc) (9) 



y't = -sin0(a;i - a;^) + cos (/)(t/i - y^). (10) 



Then the semi-major length a is determined such that 
a certain percentage. Pa, of the data points satisfy x'^ E 
[xc — a,Xc + a]. The semi-minor length b is determined in a 
similar manner with percentage Pi,. Naturally, Pa and Pi, are 
related to the noise level. If we make additional assumptions 
that the data points are independent of each other and have 
a multivariate normal distribution: z,; ^ JV{z^, A^,-), and that 
they are distributed uniformly in angle around the periphery, 
then the relationship between Pa, Pb and the noise level cr can 
be approximated as 
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and 



n(7') 



erf(-^y(l 
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de, (12) 



where 7 = - and j' — -. Recall that we have the 
approximation (j?]). Assuming that the estimated ellipse is close 
to the true model, the noise level can be readily estimated as 



1 



^ 2(1 + cos Ci)' 



(13) 



With this noise estimate, we can perform the axial guided 
fitting as described, which results in the following procedure: 

1) Solve (|5]l by gradient descent algorithms to obtain (f), 
{xc.yc), /z, and Ci, Vi; 

2) Translate and rotate the data set so that the estimated 
ellipse is located at the standard position; 

3) Estimate the noise level by ( [T3| ); 

4) Calculate Pa and Ph from (|ll|i and ( [T2] i; 

5) Find a and b so that a portion Pa of the data points 
satisfy x[ E [xc ~a,Xc + a] and a portion Pf, of the data 
points satisfy y,- E [yc^ b, yc + b]. 

Axial guided ellipse fitting divides the ellipse fitting prob- 
lem into two stages: orientation estimation and size estimation. 
In applications where the noise level is known a priori, axial 
guided fitting could be simplified and becomes more efficient 
and accurate. 

C. Weighted Objective Function 

In order to take advantage of the robust ellipse orientation 
estimation and obtain an accurate size estimation as well, we 
propose the following weighted objective function: 
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min — > — 
Ci,C2,a n ^-^ 1 



1 



- /3 cos Ci 



|x,-Ci|| + ||x,-C2||-2a)2, (14) 



ellipse is in the standard position (centered at the origin with where /3 is a tuning parameter which varies from to 1. 



When /3 = 0, (14 1 is the same as the original objective 
function, so that we can obtain accurate ellipse orientation 
estimation. On the other hand, when /? = 1, the angle 
dependent weight 1 + cos Q is applied to mimic the objective 
function of OLSF so as to obtain an accurate size estimation. 
By varying f3 from to 1 gradually, we avoid stray local 
minima at first, and aim for accuracy in the end. 

The weighted objective minimization problem is solved by 
gradient descent in such a way that l + /3cosCi, Vi is assumed 
to be constant in each iteration and updated afterwards. The 
parameter /3 is updated in a linear manner We will show the 



Original Objective Function 



Penalized Objective Function 



efficacy of this algorithm in Section VI 



V. Spheroid Fitting 

So far, we have proposed three modified algorithms based 
on the geometric definition of an ellipse. In this section, we 
generalize our method to the three-dimensional case. 

Unfortunately, general ellipsoids do not have a natural 
geometric definition similar to ellipses. Nonetheless, we can 
still generalize our algorithm to three-dimensions in the case 
of a spheroid. A spheroid is defined as a quadric surface 
obtained by rotating an ellipse about one of its principal axes; 
in other words, a spheroid is an ellipsoid with two equal semi- 
diameters. Although a spheroid is a special case of an ellipsoid, 
this special case may have broad applications. 

According to the definition, a spheroid has the same basic 
geometric property as an ellipse. This means that all algo- 
rithms proposed in Section|lIl]and|lV]can be easily generalized 
to three dimensions in the natural way. For example, in the 
case of the weighted objective function, we still have the 
optimization problem 
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/3 cos Ci 



|x,-Ci|| + ||x,-C2||-2a)2, (15) 



with Zi, Ci, C2 e M"^. 

We now have a group of algorithms almost the same as 
in the ellipse fitting case, and we will demonstrate the fitting 



result of (15 1 as an example at the end of the next section. 



VI. Experimental Results 

To demonstrate the efficacy of the algorithms described 
above, we describe a series of experiments in different settings. 
Synthetic data has been used for the simulations. 

A set of points on the perimeter of the ellipse were drawn 
according to a uniform distribution in angle. The data points 
were obtained by corrupting the true points with independent 
and identically distributed (i.i.d.) additive Gaussian noise, with 
mean and covariance cr^I. The error rate is defined as the 
normalized area difference between the true ellipse Et and the 
fitted ellipse Ef. 



where S 



error rate 



^E, U Ef ^E, fl Ef 



2SEt 



(16) 



is the area difference and 
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Fig. 2. Accuracy of the algorithms: average error rate under a vi^ide range 
of noise levels for the four proposed algorithms. The lower bound and upper 
bound of the error bars are the 20% and 80% quantiles of 50 trials. 



Convergence of the Weighted Objective Function 



Fig. 3. Convergence rate of vv'eighted objective fitting (log scale). 



A. 



denotes the area of the true ellipse. 



Comparison of the proposed algorithms 

The accuracy of the method based on the original objective 
function, penalized fitting, axial guided fitting and weighted 
objective fitting were tested and compared under a wide range 
of noise levels. 

We consider an ellipse in standard position (centered at the 
origin and without rotation) with semi-major and semi-minor 
lengths 5 and 3 respectively. Fifty data points were drawn for 
each trial for a total run of 50 trials per noise level. Fig. |2] 
shows the mean error rate for each algorithm under a range 
of noise levels (cr^ from to 0.5). The lower bounds of the 
error bars are the 20 percent quantiles and the upper bounds 
are the 80 percent quantiles of the 50 trials. 

As shown by the simulation results, the three revised 
algorithms exhibit substantial improvement compared with the 
method using the original objective function. Penalized fitting 
and weighted objective fitting have comparable performance, 
which is slightly better than that of axial guided fitting. 

Fig. |3] shows the convergence rate of weighted objective 
fitting. Note that the plot is on a log scale. So the algorithm 
converges faster than exponential in the first few iterations. As 
for penalized fitting and axial guided fitting, they converge at 
similar speeds, except that penalized fitting has an initialization 
procedure and axial guided fitting needs noise estimation. 

B. Comparison with Algebraic Fitting, OLSF, and ML 

In this subsection, we describe simulations to compare our 
method with the previous methods (algebraic fitting, OLSF 
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Fig. 4. Average error rate under a wide range of noise levels for algebraic 
fitting, OLSF, FNS and penalized fitting 



and ML), in terms of accuracy, under a high-eccentricity 
situation over a wide range of noise levels. Penalized fitting 
is used as a representative for the proposed algorithms. For 
algebraic fitting, OLSF and ML algorithms, we implemented 
the numerically stable version of algebraic fitting from [8], 
the OLSF algorithm proposed by Ahn et al. from [2], and the 
FNS ellipse fitting algorithm from [5] respectively. 

An ellipse in standard position with semi-major and semi- 
minor lengths 8 and 2, respectively, were used in the simula- 
tions. Fifty data points were drawn for each trial for a total 
run of 50 trials per noise level for penalized fitting, algebraic 
fitting and OLSF, and a total run of 1000 trials for FNS. 

Fig. |4] shows the mean error rate under a range of noise 
levels (a^ from to 0.8) for the four algorithms. Although 
penalized fitting performs slightly worse than the other three 
algorithms when the noise level is low, it outperforms them 
when the noise level increases, which shows the robustness 
of our algorithms. The curve with triangle markers represents 
the FNS algorithm. It has high accuracy when the noise level 
is low, yet it breaks down when there are data points near 
the center of the estimated ellipse, which happens often for 
moderate or high noise levels. The FNS curve represrents the 
average error on those trials (out of 1000 trials) for which 
the algorithm produced an estimate. The dotted segment of 
the FNS plot indicates that the algorithm failed to produce an 
estimate for more than 90% of the trials. 

Fig. |5]shows the comparison with error bars for the algebraic 
fitting, OLSF and penalized fitting. As in the previous case, 
the lower bounds and upper bounds of the error bars are the 
20% and 80% quantiles of 50 trials. This demonstrates the 
robustness of our algorithms. 

As for the computational cost, our algorithms perform 
almost the same as the ML algorithms and are much more 
efficient than the OLSF algorithms in a typical setting. 

C. Spheroid Fitting 

Here we present an example of a typical spheroid fitting 
result to demonstrate the efficacy of our spheroid fitting 
algorithm based on the weighted objective function. Fifty 
true points are generated from the surface ofal0x2x2 
spheroid centered at the origin with rotation. The data points 
are generated from the true points with additive Gaussian noise 




Fig. 5. Average error rate for algebraic fitting. OLSF, penalized fitting with 
20% quantile and 80% quantile error bars. 



Spheroid Fitting (N = 100. a= 10. b = c = 2. = 0.2) 




Fig. 6. Spheroid fitting result 

with mean and variance 0.21. Fig. [6] shows the fitting result 
of our algorithm. 
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